# Clean Environment -------------------------------------------------------
rm(list = ls())
source("figure_design.R")

# TABLE B.1 ---------------------------------------------------------------

# Load the measurement model data
load("measurement_model.RData")

# Print the parameter estimates to a LaTeX table
print(xtable(parameterEstimates(fit_ord)), file = "tab_b1.tex")

# FIGURE B.1 --------------------------------------------------------------

# Read the data
covariates <- read_rds("df_full_uod.rds")

# Prepare the data for the density plot
filtered_data <- covariates %>%
  dplyr::select(lib, maj, auth) %>%
  pivot_longer(cols = everything(), names_to = "variable", values_to = "value") %>%
  mutate(variable = dplyr::recode(variable,
                           lib = "Liberal",
                           maj = "Majoritarian",
                           auth = "Authoritarian"))

# Create the density plot
density_plot <- ggplot(filtered_data, aes(x = value, fill = variable)) +
  geom_density(alpha = 0.5) +
  facet_wrap(~ variable, scales = "free") +
  xlab("") +
  ylab("Density") +
  scale_fill_manual(values = plot_colors) +
  theme(legend.position = "none")

# Save the plot to a PDF file
pdf(file = "fig_b1.pdf", width = 6, height = 2)
print(density_plot)
dev.off()

# FIGURE B.2 --------------------------------------------------------------

# Prepare the data for the density plot with party preference
filtered_data <- covariates %>%
  dplyr::select(lib, maj, auth, party_preference) %>%
  pivot_longer(cols = c(lib, maj, auth), names_to = "variable", values_to = "value") %>%
  mutate(variable = dplyr::recode(variable,
                           lib = "Liberal",
                           maj = "Majoritarian",
                           auth = "Authoritarian"))

# Calculate mean values
mean_values <- filtered_data %>%
  group_by(variable, party_preference) %>%
  summarise(mean_value = mean(value, na.rm = TRUE), .groups = 'drop')

# Create the density plot with mean values
density_plot <- ggplot(filtered_data, aes(x = value, fill = variable)) +
  geom_density(alpha = 0.5) +
  facet_grid(party_preference ~ variable, scales = "free") +
  xlab("Value") +
  ylab("Density") +
  scale_fill_manual(values = plot_colors) +
  theme(legend.position = "none") +
  geom_text(data = mean_values, aes(x = mean_value, y = 0, label = paste("Mean:", round(mean_value, 2))),
            vjust = -1.5, hjust = -0.1, size = 3, color = "black", inherit.aes = FALSE)

# Save the plot to a PDF file
pdf(file = "fig_b2.pdf", width = 8, height = 12)
print(density_plot)
dev.off()

# FIGURE B.3 --------------------------------------------------------------

# Read the data
df_uod <- read_rds("df_uod_items.rds") %>% dplyr::select(-ResponseId) %>%
  mutate(across(everything(), as.numeric))

# Perform parallel analysis to retrieve the optimal number of factors
fa.parallel(df_uod, fa = "fa")

# Save the parallel analysis plot
dev.copy(pdf, "fig_b3.pdf", width = 12, height = 6)
dev.off()

# TABLE B.2 ---------------------------------------------------------------

# Perform EFA with a 3-factor solution using Promax rotation
efa <- factanal(df_uod, 3, rotation = "promax")
print(efa, digits = 2, cutoff = .2, sort = TRUE)

# Extract loadings and uniquenesses
loadings <- as.data.frame(efa$loadings[, 1:3])
uniquenesses <- efa$uniquenesses

# Combine the results into a single data frame
results <- cbind(loadings, uniquenesses)
colnames(results) <- c("Factor1", "Factor2", "Factor3", "Uniquenesses")

# Convert the results to a LaTeX table
results_xtable <- xtable(results, caption = "Factor Analysis Results (Promax Rotation)", label = "tab:factor_analysis_promax")

# Print the LaTeX code
print(results_xtable, type = "latex", include.rownames = TRUE, booktabs = TRUE, file = "tab_b2.tex")

# TABLE B.4 ---------------------------------------------------------------

# Perform EFA with a 3-factor solution using Varimax rotation
efa_varimax <- factanal(df_uod, 3, rotation = "varimax")
print(efa_varimax, digits = 2, cutoff = .3, sort = TRUE)

# Extract loadings and uniquenesses
loadings_varimax <- as.data.frame(efa_varimax$loadings[, 1:3])
uniquenesses_varimax <- efa_varimax$uniquenesses

# Combine the results into a single data frame
results_varimax <- cbind(loadings_varimax, uniquenesses_varimax)
colnames(results_varimax) <- c("Factor1", "Factor2", "Factor3", "Uniquenesses")

# Convert the results to a LaTeX table
results_xtable_varimax <- xtable(results_varimax, caption = "Factor Analysis Results (Varimax Rotation)", label = "tab:factor_analysis_varimax")

# Print the LaTeX code
print(results_xtable_varimax, type = "latex", include.rownames = TRUE, booktabs = TRUE, file = "tab_b3.tex")